Εργαστηριακή Άσκηση 1#

Σκοπός της πρώτης σειράς ασκήσεων είναι, αφ’ ενός η εξοικείωση με το προγραμματιστικό περιβάλλον της Python, αφ’ ετέρου, η εισαγωγή στους τρόπους παράστασης και επεξεργασίας τηλεπικοινωνιακών σημάτων στη συγκεκριμένη γλώσσα προγραμματισμού.

Μέρος 1: Εξοικείωση με το προγραμματιστικό περιβάλλον της Python#

Λίγα λόγια για τον interactive python interpreter [ipython] και το jupyter notebook

Interactive python interpreter (IPython) notebook#

Το IPython (Interactive Python) (https://ipython.org/) είναι ένα διαδραστικό κέλυφος γραμμής εντολών σχεδιασμένο για τη γλώσσα προγραμματισμού Python.Υπερβαίνει τον παραδοσιακό διερμηνέα Python και προσφέρει βελτιωμένες δυνατότητες (REPL):

1)Ανάγνωσης (Read)
2)Aξιολόγησης (Evaluate)
3)Eκτύπωσης (Print)
4)Επαναλάψης (Loop)

Jupyter notebook#

To Project Jupyter είναι ένας μη κερδοσκοπικός οργανισμός με αποστολή τη συγγραφή ανοικτού λογισμικού για διαδραστικές εφαρμογές. Ξεκίνησε από ipython, ωστόσο σήμερα προσφέρει προγράμματα σε πολλές γλώσσες προγραμματισμού.

Το jupyter notebook (https://jupyter.org/) είναι μια πλατφόρμα web ανοικτού λογισμικού για την ανάπτυξη διαδραστικών εφαρμογών, κυρίως για επεξεργασία (επιστημονικών) δεδομένων και μηχανικής μάθησης.

Εξάσκηση#

from scipy import signal
# Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal
# https://docs.scipy.org/doc/scipy/reference/signal.html
from scipy.fft import fft, fftfreq
from dash import Dash, dcc, html, Input, Output
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
import pandas
import plotly
import plotly.express as px
import plotly.graph_objs as go
# Δείτε την έκδοση της αριθμητικής βιβλιοθήκης numpy
import dash_bootstrap_components as dbc
np.__version__
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 5
      2 # Ανατρέξτε στην τεκμηρίωση της βιβλιοθήκης scipy.signal
      3 # https://docs.scipy.org/doc/scipy/reference/signal.html
      4 from scipy.fft import fft, fftfreq
----> 5 from dash import Dash, dcc, html, Input, Output
      6 import matplotlib
      7 import matplotlib.pyplot as plt

ModuleNotFoundError: No module named 'dash'
import warnings
warnings.filterwarnings('ignore')

Μην ξεχνάτε ότι η IPython μας δίνει τη δυνατότητα να ‘εξερευνήσουμε’ το περιεχόμενο ενός package, χρησιμοποιώντας τη δυνατότητα του tab-completion, ή τη χρήση του ? για help/documentation: Π.χ., για να δούμε όλα τα περιεχόμενα του signal namespace δίνουμε:

In [3]: signal?

και για να καλέσουμε την ενσωμετωμένη τεκμηρίωση της numpy, δίνουμε:

In [4]: np?

Περισσότερες πληροφορίες μπορείτε να πάρετε από το http://www.numpy.org.

# Δημιουργήστε ένα βαθμωτό (μονοδιάστατο) μέγεθος

s=2
print('s=',s)
s= 2
# Δημιουργείστε ένα διάνυσμα πραγματικών τιμών:
# Στο MATLAB: v = [1,5,9] ή v = [1 5 9]

v=np.array([1,5,9])
print('v=',v)
v= [1 5 9]
# Δημιουργείστε έναν πίνακα πραγματικών τιμών:
# Στο MATLAB: a = a=[[1,2,3];[4,5,6];[7,8,9]] ή a=[1,2,3;4,5,6;7,8,9]
a=np.array([[1,2,3],[4,5,6],[7,8,9]])
print('a=',a)
a= [[1 2 3]
 [4 5 6]
 [7 8 9]]
# Αθροίστε

a+5
array([[ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])
#Πολλαπλασιάστε

b=s*v*2
print('b=',b)
b= [ 4 20 36]
# Πολλαπλασιάστε στοιχείο-προς-στοιχείο (elementwise)
# MATLAB: v.*b

np.multiply(v,b)
array([  4, 100, 324])
#Ελέγξτε το μήκος ενός διανύσματος
# MATLAB: length(v)

len(v)
3
# Ελέγξτε το μέγεθος ενός πίνακα
# MATLAB: size(a)

a.shape   # για array: np.array(a.shape)
(3, 3)
# Προσπελάστε συγκεκριμένα στοιχεία ενός πίνακα
# Η δεικτοδότηση αρχίζει από το 0. 
# MATLAB: a(1,2)  
# --- ΠΡΟΣΟΧΗ, στο MATLAB η δεικτοδ΄ότηση αρχ΄ίζει από το 1!

a[0,1]
2
# Προσπελάστε συγκεκριμένα στοιχεία ενός πίνακα (συνέχεια)
# Αρνητικές τιμές μετρούν από το τέλος, π.χ. το -1 
# αναφέρεται στο τελευταίο στοιχείο)

a[1,-1]
6
# Προσπελάστε συγκεκριμένο τμήμα ενός διανύσματος
# MATLAB: v(1:9)

v[1:3]

# ΠΡΟΣΟΧΗ: τα στοιχεία [2ο,3ο] δίνονται ως 1:3 και όχι ως 1:2
# Δοκιμάστε το v[1:2]...
array([5, 9])
# Προσπελάστε συγκεκριμένα τμήματα ενός πίνακα

a[0:2,:]

# Ομοίως: οι γραμμές 1 & 2 δίνονται ως 0:2 και όχι ως 0:1
array([[1, 2, 3],
       [4, 5, 6]])
# Δημιουργήστε ένα διάνυσμα με στοιχεία από το 0 έως το 0.5 και βήμα 0.1
# MATLAB: t=(0:0.1:0.4)

t=np.arange(0,0.5,0.1)
print('t=',t)
t= [0.  0.1 0.2 0.3 0.4]

Μέρος 2: Δειγματοληψία - Ψηφιοποίηση#

Τα πρωτογενή σήματα είναι κυρίως αναλογικά (συνεχούς χρόνου). Για να τα παραστήσουμε και επεξεργαστούμε στον υπολογιστή μας (ή άλλη ψηφιακή μηχανή) θα πρέπει πρώτα να τα ψηφιοποιήσουμε. Υποθέστε ένα σήμα συνεχούς χρόνου \(x(t)\) με μετασχηματισμό Fourier (Continuous Time Fourier Transform – CTFT): $\(X(f)=\int_{-\infty}^{\infty} x(t)e^{-j2\pi ft} dt\)$

app1 = Dash(__name__, external_stylesheets=[dbc.themes.BOOTSTRAP])

# Παράμετροι δειγματοληψίας
fs = 1000  # Sampling frequency in Hz
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

# Layout of the Dash app
app1.layout = html.Div(
    children=[
        html.H2('Interactive Signal and Fourier Transformation'),
        dcc.RangeSlider(
            id='frequency-slider',
            min=1,
            max=51,
            step=1,
            value=[5, 10],
            marks={i: f'{i} Hz' for i in range(1, 51, 5)}
        ),
        dbc.Row(
            children=[
                dbc.Col(
                    dcc.Graph(id='original-signal'),
                    width=6
                ),
                dbc.Col(
                    dcc.Graph(id='fourier-transform'),
                    width=6
                )
            ]
        )
    ],
    style={'backgroundColor': 'white', 'padding': '20px'}
)

# Callback to update graphs
@app1.callback(
    [Output('original-signal', 'figure'),
     Output('fourier-transform', 'figure')],
    [Input('frequency-slider', 'value')]
)
def update_graph(selected_frequencies):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(selected_frequencies[0], selected_frequencies[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
                    

    # Calculate the FFT of the original signal
    original_signal_fft = fft(signal)
    # Calculate frequencies for the FFT of the original signal
    original_freqs = fftfreq(L, T)[:L // 2]
    # Calculate magnitude of Fourier coefficients (amplitude) for the original signal
    original_magnitude = np.abs(original_signal_fft)[:L // 2]
   

    # Original signal graph
    original_signal_figure = {
        'data': [go.Scatter(x=t, y=signal, mode='lines', line=dict(color='#00CC96', width=2))],
        'layout': go.Layout(
            title='Original Signal',
            xaxis={'title': 'Time (s)', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Fourier Transformation graph
    fourier_transform_figure = {
        'data': [go.Scatter(x=original_freqs, y=original_magnitude, mode='lines', line=dict(color='#1F77B4', width=2))],
        'layout': go.Layout(
            title='Fourier Transformation of Original Signal',
            xaxis={'title': 'Frequency (Hz)', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    
    # Return this additional figure in the callback
    return original_signal_figure, fourier_transform_figure

    

if __name__ == '__main__':
    app1.run_server(debug=True, port=8050)

Λαμβάνοντας δείγματα του \(x(t)\) με ρυθμό \(f_s=1/T_s\) παράγεται σήμα διακριτού χρόνου \(x(nT_s)\). Μαθηματικά το αναπαριστάνουμε ως σειρά συναρτήσεων δέλτα $\(x_\delta (t)=\sum_{n=-\infty}^{\infty}x(nT_s)\delta(t-nT_s)=x(t)\sum_{n=-\infty}^{\infty}\delta(t-nT_s)\)\( με μετασχηματισμό Fourier \)\(X_\delta (f)=\sum_{n=-\infty}^{\infty}x(nT_s)e^{-j2\pi fnT_s}=X(f)*1/T_s\sum_{n=-\infty}^{\infty}\delta(f-k/T_s)=1/T_s\sum_{n=-\infty}^{\infty}X(f-k/T_s)\)$ που είναι περιοδική συνάρτηση.

from scipy.fftpack import fftshift, ifftshift


app2 = Dash(__name__, external_stylesheets=[dbc.themes.BOOTSTRAP])

# Παράμετροι δειγματοληψίας
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

# Layout of the Dash app
app2.layout = html.Div(
    children=[
        html.H2('Interactive Signal and Fourier Transformation'),
        dcc.Slider(
            id='samples-slider',
            min=100,
            max=1000,
            step=25,
            value=500, 
            marks={i: f'{i} samples' for i in range(100, 1000, 100)},
            included=False,  # To remove the color fill
            updatemode='drag',  # Update slider continuously while dragging
            tooltip={'placement': 'bottom'}  # Show tooltip at the bottom
        ),
        dcc.RangeSlider(
            id='frequency-slider',
            min=1,
            max=51,
            step=1,
            value=[5, 10],
            marks={i: f'{i} Hz' for i in range(1, 51, 5)}
        ),
        dbc.Row(
            children=[
                dbc.Col(
                    dcc.Graph(id='sampled-signal'),
                    width=6
                ),
                dbc.Col(
                    dcc.Graph(id='sampled-fourier-signal'),
                    width=6
                )

            ]
        )
    ],
    style={'backgroundColor': 'white', 'padding': '20px'}
)

# Callback to update graphs
@app2.callback(
    [Output('sampled-signal', 'figure'),
     Output('sampled-fourier-signal', 'figure')],
    [Input('frequency-slider', 'value'),
    Input('samples-slider', 'value')]
)
def update_graph(selected_frequencies,selected_samples):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(selected_frequencies[0], selected_frequencies[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
    
   # Correct downsampling approach
    sample_rate = selected_samples  # Directly using selected_samples as an integer
    downsampling_factor = int(fs / sample_rate)  # Downsampling factor

    # Select every nth sample from the original time vector to match the downsampled signal
    sample_points = t[::downsampling_factor]

    # Downsample the signal by selecting every nth sample
    downsampled_signal = signal[::downsampling_factor]

    # Recalculate L for the downsampled signal if needed
    L_downsampled = len(downsampled_signal)

    # Calculate the DFT of the downsampled signal
    sampled_signal_fft = fft(downsampled_signal)

    # Calculate frequencies for the FFT of the downsampled signal
    # Note: The new sampling period is the inverse of the new sampling rate
    sampled_freqs = fftfreq(L_downsampled, 1/sample_rate)[:L_downsampled // 2]

    # Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signal
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    n = len(downsampled_signal)
    T = 1 / sample_rate  # Recalculate the sampling period with the selected sample rate

    # After calculating the magnitude of Fourier coefficients
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    # Create a new frequency vector that includes the replicated frequencies
    # First, calculate the original frequency bins for the positive frequencies
    original_sampled_freqs = fftfreq(L_downsampled, T)[:L_downsampled // 2]

    shifted_freqs_negative = original_sampled_freqs - 1/T  # Shift for -1/Ts
    shifted_freqs_positive = original_sampled_freqs + 1/T  # Shift for +1/Ts

    # Concatenate the original and shifted (replicated) frequencies and magnitudes
    # This includes the original frequencies, and the replications at -1/Ts and +1/Ts
    replicated_freqs = np.concatenate([shifted_freqs_negative, original_sampled_freqs, shifted_freqs_positive])
    replicated_magnitude = np.concatenate([sampled_magnitude, sampled_magnitude, sampled_magnitude])

    # Sort the replicated frequencies and magnitudes in ascending order for proper plotting
    sorted_indices = np.argsort(replicated_freqs)
    sorted_replicated_freqs = replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Sampled signal graph
    sampled_signal_figure = {
        'data': [go.Scatter(x=sample_points, y=downsampled_signal, mode='markers', marker=dict(color='#00CC96', size=7))],
        'layout': go.Layout(
            title='Sampled Signal',
            xaxis={'title': 'Time (s)', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Fourier Transformation graph for the sampled signal
    sampled_fourier_transform_figure = {
        'data': [go.Scatter(x=sorted_replicated_freqs, y=sorted_replicated_magnitude, mode='lines', line=dict(color='#1F77B4', width=2))],
        'layout': go.Layout(
            title='Fourier Transformation of Sampled Signal',
            xaxis={'title': 'Frequency (Hz)', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Return this additional figure in the callback
    return sampled_signal_figure, sampled_fourier_transform_figure

    
    

if __name__ == '__main__':
    app2.run_server(debug=True, port=8051)

Για βαθυπερατά σήματα \(x(t)\) εύρους ζώνης W, με την υπόθεση ότι ο ρυθμός δειγματοληψίας \(fs ≥ 2W\), ισχύει ότι \(X(f) = T_s X_\delta(f)\), \(0 ≤ f ≤ W\), δηλαδή, το σήμα \(X(f)\) προκύπτει μετά από διάβαση του δειγματοληπτημένου \(x_\delta(t)\) μέσω ιδανικού βαθυπερατού φίλτρου κέρδους \(T_s\). Από το προηγούμενο σχήμα γίνεται φανερό ότι εάν η δειγματοληψία γίνει με συχνότητα μικρότερη του διπλασίου της ανώτερης συχνότητας \(W\) του σήματος (υποδειγμάτιση – undersampling), τότε εμφανίζονται στην περιοχή συχνοτήτων του σήματος «είδωλα» φάσματος από ανώτερες συχνότητες που δεν επιτρέπουν την ακριβή αποκατάσταση του αρχικού σήματος συνεχούς χρόνου. Το φαινόμενο αυτό ονομάζεται αναδίπλωση ή επικάλυψη (aliasing), το δε σφάλμα κατά την αποκατάσταση του αρχικού σήματος αποκαλείται σφάλμα αναδίπλωσης (aliasing error). Η δειγματοληψία στο πεδίο του χρόνου αποτελεί τη βάση για τον ορισμό του μετασχηματισμού Fourier διακριτού χρόνου (Discrete Time Fourier Transform – DTFT). Για μια σειρά διακριτών αριθμών \(x[n]\), ο μετασχηματισμός Fourier διακριτού χρόνου ορίζεται ως:

3.PNG

O DTFT είναι περιοδική συνάρτηση με περίοδο \(1\), επομένως, αρκεί ο υπολογισμός του στο διάστημα συχνοτήτων \([0,1]\) ή ισοδύναμα \([-½,½]\). Να σημειωθεί ότι ο DTFT, παρότι προκύπτει από μια σειρά διακριτών αριθμών \(x[n]\), είναι συνεχής συνάρτηση της μεταβλητής \(\phi\) όπως παραστατικά φαίνεται στο επόμενο σχήμα.

from scipy.fftpack import fftshift, ifftshift


app3 = Dash(__name__, external_stylesheets=[dbc.themes.BOOTSTRAP])

# Παράμετροι δειγματοληψίας
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

# Layout of the Dash app
app3.layout = html.Div(
    children=[
        html.H2('Interactive Signal and Fourier Transformation'),
        dcc.Slider(
            id='samples-slider',
            min=100,
            max=1000,
            step=25,
            value=500, 
            marks={i: f'{i} samples' for i in range(100, 1000, 100)},
            included=False,  # To remove the color fill
            updatemode='drag',  # Update slider continuously while dragging
            tooltip={'placement': 'bottom'}  # Show tooltip at the bottom
        ),
        dcc.RangeSlider(
            id='frequency-slider',
            min=1,
            max=51,
            step=1,
            value=[5, 10],
            marks={i: f'{i} Hz' for i in range(1, 51, 5)}
        ),
        dbc.Row(
            children=[
                dbc.Col(
                    dcc.Graph(id='sampled-signal'),
                    width=6
                ),
                dbc.Col(
                    dcc.Graph(id='sampled-fourier-signal'),
                    width=6
                )

            ]
        )
    ],
    style={'backgroundColor': 'white', 'padding': '20px'}
)

# Callback to update graphs
@app3.callback(
    [Output('sampled-signal', 'figure'),
     Output('sampled-fourier-signal', 'figure')],
    [Input('frequency-slider', 'value'),
    Input('samples-slider', 'value')]
)
def update_graph(selected_frequencies,selected_samples):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(selected_frequencies[0], selected_frequencies[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
    
   # Correct downsampling approach
    sample_rate = selected_samples  # Directly using selected_samples as an integer
    downsampling_factor = int(fs / sample_rate)  # Downsampling factor

    # Downsample the signal by selecting every nth sample
    downsampled_signal = signal[::downsampling_factor]

    # Recalculate L for the downsampled signal if needed
    L_downsampled = len(downsampled_signal)

    # Calculate the DFT of the downsampled signal
    sampled_signal_fft = fft(downsampled_signal)

    # Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signal
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    n = len(downsampled_signal)
    T = 1 / sample_rate  # Recalculate the sampling period with the selected sample rate

    # After calculating the magnitude of Fourier coefficients
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    # Create a new frequency vector that includes the replicated frequencies
    # First, calculate the original frequency bins for the positive frequencies
    original_sampled_freqs = fftfreq(L_downsampled, T)[:L_downsampled // 2]

    shifted_freqs_negative = original_sampled_freqs - 1/T  # Shift for -1/Ts
    shifted_freqs_positive = original_sampled_freqs + 1/T  # Shift for +1/Ts

    # Concatenate the original and shifted (replicated) frequencies and magnitudes
    # This includes the original frequencies, and the replications at -1/Ts and +1/Ts
    replicated_freqs = np.concatenate([shifted_freqs_negative, original_sampled_freqs, shifted_freqs_positive])
    replicated_magnitude = np.concatenate([sampled_magnitude, sampled_magnitude, sampled_magnitude])

    # Sort the replicated frequencies and magnitudes in ascending order for proper plotting
    sorted_indices = np.argsort(replicated_freqs)
    sorted_replicated_freqs = replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Generate an array of sample indices for the downsampled signal
    sample_indices = np.arange(len(downsampled_signal))

    # Update the Sampled signal graph to plot against sample indices
    sampled_signal_figure = {
        'data': [go.Scatter(x=sample_indices, y=downsampled_signal, mode='markers', marker=dict(color='#00CC96', size=7))],
        'layout': go.Layout(
            title='Sampled Signal',
            xaxis={'title': 'Sample Number', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Normalize the frequency values
    normalized_freqs = original_sampled_freqs / (1/T)

    # Since you previously concatenated and sorted for replication, ensure to apply normalization there as well
    normalized_replicated_freqs = np.concatenate([
        (shifted_freqs_negative / (1/T)),  # Normalize -1/Ts shifted frequencies
        normalized_freqs,  # Already normalized original frequencies
        (shifted_freqs_positive / (1/T))   # Normalize +1/Ts shifted frequencies
    ])

    # Sort the normalized and replicated frequencies for proper plotting
    sorted_indices = np.argsort(normalized_replicated_freqs)
    sorted_normalized_replicated_freqs = normalized_replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Update the Fourier Transformation graph with normalized frequencies
    sampled_fourier_transform_figure = {
        'data': [go.Scatter(x=sorted_normalized_replicated_freqs, y=sorted_replicated_magnitude, mode='lines', line=dict(color='#1F77B4', width=2))],
        'layout': go.Layout(
            title='Normalized Fourier Transformation of Sampled Signal',
            xaxis={'title': 'Normalized Frequency (f/fs)', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Return this additional figure in the callback
    return sampled_signal_figure, sampled_fourier_transform_figure

    
    

if __name__ == '__main__':
    app3.run_server(debug=True, port=8051)

Με τη σειρά των διακριτών αριθμών να προκύπτει ως αποτέλεσμα δειγματοληψίας, \(x[n]=x(nT_s)\), ο DTFT και ο μετασχηματισμός Fourier \(X_\delta(f)\) του δειγματοληπτημένου σήματος συνδέονται μέσω της αντιστοιχίας \(\phi ↔ f/f_s\). Η συνήθης πρακτική είναι να παριστάνουμε τον λόγο \(f/f_s\) ως κανονικοποιημένη συχνότητα \(\phi\) (\(f_D\), στις σημειώσεις σας) και οι πραγματικές συχνότητες να προκύπτουν ως πολλαπλάσιά της (συνήθως κλασματικά). Για τη σύνδεση του DTFT με τον μετασχηματισμό Fourier \(X(f)\) του σήματος πρέπει επιπλέον να γίνει αναγωγή στην περίοδο δειγματοληψίας με πολλαπλασιασμό επί \(T_s\) (ή διαίρεση με \(f_s\)). Κατ΄ αναλογία με τη δειγματοληψία σημάτων στο χρόνο μπορούμε να κάνουμε δειγματοληψία στο πεδίο της συχνότητας λαμβάνοντας διακριτές τιμές \(X(kf_o)\) του μετασχηματισμού Fourier που αντιστοιχούν σε ανάλυση συχνότητας \(f_o=1/T_o\). Αυτό ισοδυναμεί με περιοδική επανάληψη του σήματος συνεχούς χρόνου \(x(t)\) κάθε \(Τ_ο\), αφού το περιοδικό σήμα

5.PNG

έχει μετασχηματισμό Fourier

6.PNG

Επομένως, \(X[k] = X(kf_o)/Τ_o\) είναι οι συντελεστές του αναπτύγματος σε σειρά Fourier.του περιοδικού σήματος \(x_p(t)\). Προφανώς, για σήματα \(x(t)\) πεπερασμένης διάρκειας, όπου \(x(t)=0\) για \(|t| ≥ T\), με την υπόθεση ότι η περίοδος \(T_o ≥ 2T\), ισχύει ότι \(x(t) = x_p(t)\) για \(|t| ≤ T\). Στην πράξη, τα σήματα έχουν πολύ μεγάλη διάρκεια για να μπορέσουμε να τα αναλύσουμε στην ολότητά τους. Έτσι εφαρμόζουμε ένα ορθογωνικό χρονικό παράθυρο, ώστε να διατηρήσουμε μόνο το πιο σημαντικό τους μέρος για το διάστημα παρατήρησης και \(x(t)= 0\), αλλού. Κατά τον υπολογισμό του DTFT \(X_d(\phi)\) ενός τέτοιου ακρωτηριασμένου σήματος, αντί του απείρου αθροίσματος, περιοριζόμαστε σε μια πεπερασμένου μήκους \(L\) σειρά αριθμών \(x[n]\), οπότε

1.PNG

H δειγματοληψία του \(X_d(\phi)\) στο πεδίο συχνότητας σε \(Ν\) ισαπέχουσες κανονικοποιημένες συχνότητες \(0\), \(1/Ν\), \(2/Ν\), \(…\), \((Ν-1)/Ν\), δίνει

1.PNG

όπου, εάν \(N≥L\), θέτουμε \(x[n]=0\) για \(n≥L\). Η τελευταία σχέση αναγνωρίζεται ως ο διακριτός μετασχηματισμός Fourier (Discrete Fourier Transform – DFT), ο οποίος για μια πεπερασμένη σειρά \(xn\), \(n=0\), \(1\), \(…\), \(N-1\), ορίζεται ως:

1.PNG

και ο αντίστροφός του είναι

1.PNG

Η \(X_d(\phi)\) ως DTFT είναι περιοδική συνάρτηση και εάν η αρχική σειρά xn ήταν περιοδική (και δεν εφαρμόζαμε το παράθυρο), τότε η \(X_d(\phi)\) θα ήταν μηδέν παντού εκτός των σημείων της δειγματοληψίας \(k/Ν\). Δηλαδή, εάν θεωρήσουμε μια πεπερασμένου μήκους σειρά αριθμών που επαναλαμβάνεται περιοδικά, o διακριτού χρόνου μετασχηματισμός Fourier της (DTFT) είναι και αυτός περιοδικός και διακριτός. Επιπλέον, ο DFT και ο αντίστροφός του IDFT, εάν δεν περιορίζαμε τους δείκτες \(n\) και \(k\) μεταξύ \(0\) και \(N-1\), θα ήταν περιοδικές συναρτήσεις. Άρα η πεπερασμένη σειρά xn μπορεί να θεωρηθεί ως ένα περιοδικό σήμα διακριτού χρόνου ιδωμένο μόνο κατά τη διάρκεια μιας περιόδου και ο DFT, η σειρά \(X_k\), ως τα δείγματα με ανάλυση \(1/Ν\) του DTFT \(X_d(\phi)\) στο πεδίο κανονικοποιημένων συχνοτήτων \([0,1]\), όπως φαίνεται στο επόμενο σχήμα.

from scipy.fftpack import fftshift, ifftshift


app4 = Dash(__name__, external_stylesheets=[dbc.themes.BOOTSTRAP])

# Παράμετροι δειγματοληψίας
T = 1 / fs  # Sampling period
L = 1000  # Number of samples

# Create time vector
t = np.arange(0, L) * T

# Layout of the Dash app
app4.layout = html.Div(
    children=[
        html.H2('Interactive Signal and Fourier Transformation'),
        dcc.Slider(
            id='samples-slider',
            min=100,
            max=1000,
            step=25,
            value=500, 
            marks={i: f'{i} samples' for i in range(100, 1000, 100)},
            included=False,  # To remove the color fill
            updatemode='drag',  # Update slider continuously while dragging
            tooltip={'placement': 'bottom'}  # Show tooltip at the bottom
        ),
        dcc.RangeSlider(
            id='frequency-slider',
            min=1,
            max=51,
            step=1,
            value=[5, 40],
            marks={i: f'{i} Hz' for i in range(1, 51, 5)}
        ),
        dcc.Slider(
            id='N-slider',
            min=1,
            max=25,
            step=1,
            value=1, 
            marks={i: f'{i} N' for i in range(1, 25, 1)},
            included=False,  # To remove the color fill
            updatemode='drag',  # Update slider continuously while dragging
            tooltip={'placement': 'bottom'}  # Show tooltip at the bottom
        ),
        dbc.Row(
            children=[
                dbc.Col(
                    dcc.Graph(id='sampled-signal'),
                    width=6
                ),
                dbc.Col(
                    dcc.Graph(id='sampled-fourier-signal'),
                    width=6
                )
            ]
        )
    ],
    style={'backgroundColor': 'white', 'padding': '20px'}
)

# Callback to update graphs
@app4.callback(
    [Output('sampled-signal', 'figure'),
     Output('sampled-fourier-signal', 'figure')],
    [Input('frequency-slider', 'value'),
    Input('samples-slider', 'value'),
    Input('N-slider','value')]
)
def update_graph(selected_frequencies,selected_samples,selected_N):
    # Generate the signal with the selected frequencies
    signal = np.zeros(len(t))
    for f in range(selected_frequencies[0], selected_frequencies[1] + 1):
        signal += np.sin(2 * np.pi * f * t)
    
   # Correct downsampling approach
    sample_rate = selected_samples  # Directly using selected_samples as an integer
    downsampling_factor = int(fs / sample_rate)  # Downsampling factor

    # Downsample the signal by selecting every nth sample
    downsampled_signal = signal[::downsampling_factor]

    # Recalculate L for the downsampled signal if needed
    L_downsampled = len(downsampled_signal)

    # Calculate the DFT of the downsampled signal
    sampled_signal_fft = fft(downsampled_signal)

    # Calculate magnitude of Fourier coefficients (amplitude) for the downsampled signal
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    n = len(downsampled_signal)
    T = 1 / sample_rate  # Recalculate the sampling period with the selected sample rate

    # After calculating the magnitude of Fourier coefficients
    sampled_magnitude = np.abs(sampled_signal_fft)[:L_downsampled // 2]

    # Create a new frequency vector that includes the replicated frequencies
    # First, calculate the original frequency bins for the positive frequencies
    original_sampled_freqs = fftfreq(L_downsampled, T)[:L_downsampled // 2]

    shifted_freqs_negative = original_sampled_freqs - 1/T  # Shift for -1/Ts
    shifted_freqs_positive = original_sampled_freqs + 1/T  # Shift for +1/Ts

    # Concatenate the original and shifted (replicated) frequencies and magnitudes
    # This includes the original frequencies, and the replications at -1/Ts and +1/Ts
    replicated_freqs = np.concatenate([shifted_freqs_negative, original_sampled_freqs, shifted_freqs_positive])
    replicated_magnitude = np.concatenate([sampled_magnitude, sampled_magnitude, sampled_magnitude])

    # Sort the replicated frequencies and magnitudes in ascending order for proper plotting
    sorted_indices = np.argsort(replicated_freqs)
    sorted_replicated_freqs = replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Generate an array of sample indices for the downsampled signal
    sample_indices = np.arange(len(downsampled_signal))

    # Create data for the markers on top of the stems
    sampled_signal_markers = go.Scatter(
        x=sample_indices, 
        y=downsampled_signal, 
        mode='markers', 
        marker=dict(color='#00CC96', size=5),
        name='Markers'
    )

    # Create a list to hold all the scatter plots for the stems
    sampled_signal_stems = []

    # For each sample point, create a line from the x-axis to the point
    for x, y in zip(sample_indices, downsampled_signal):
        sampled_signal_stems.append(
            go.Scatter(
                x=[x, x],
                y=[0, y],
                mode='lines',
                line=dict(color='#00CC96', width=0.3),
                showlegend=False  # Hide legend for each stem
            )
        )

    # Combine markers and stems for the plot
    sampled_signal_data = [sampled_signal_markers] + sampled_signal_stems

    # Update the figure with the new layout
    sampled_signal_figure = {
        'data': sampled_signal_data,
        'layout': go.Layout(
            title='Sampled Signal',
            xaxis={'title': 'Sample Number', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Normalize the frequency values
    normalized_freqs = original_sampled_freqs / (1/T)

    # Since you previously concatenated and sorted for replication, ensure to apply normalization there as well
    normalized_replicated_freqs = np.concatenate([
        (shifted_freqs_negative / (1/T)),  # Normalize -1/Ts shifted frequencies
        normalized_freqs,  # Already normalized original frequencies
        (shifted_freqs_positive / (1/T))   # Normalize +1/Ts shifted frequencies
    ])

    # Sort the normalized and replicated frequencies for proper plotting
    sorted_indices = np.argsort(normalized_replicated_freqs)
    sorted_normalized_replicated_freqs = normalized_replicated_freqs[sorted_indices]
    sorted_replicated_magnitude = replicated_magnitude[sorted_indices]

    # Let's assume N is defined; for example, N could be 10 for sampling every 1/10th
    N=selected_N

    # Adjust the sampling interval for selecting frequencies and magnitudes
    sampled_indices = np.arange(0, len(sorted_normalized_replicated_freqs), N)
    sampled_normalized_replicated_freqs = sorted_normalized_replicated_freqs[sampled_indices]
    sampled_replicated_magnitude = sorted_replicated_magnitude[sampled_indices]

    markers_data = go.Scatter(
    x=sampled_normalized_replicated_freqs, 
    y=sampled_replicated_magnitude, 
    mode='markers', 
    marker=dict(color='#00CC96', size=5),
    name='Markers',
    showlegend=False
)

    # Create a list to hold all the scatter plots for the stems
    stems_data = []

    # For each point, create a line from the x-axis to the point
    for x, y in zip(sampled_normalized_replicated_freqs, sampled_replicated_magnitude):
        stems_data.append(
            go.Scatter(
                x=[x, x],
                y=[0, y],
                mode='lines',
                line=dict(color='#00CC96', width=0.25),
                showlegend=False  # Hide legend for each stem
            )
        )

    # Combine markers and stems for the plot
    data = [markers_data] + stems_data

    # Update the figure with the new layout
    sampled_fourier_transform_figure = {
        'data': data,
        'layout': go.Layout(
            title='Sampled Normalized Fourier Transformation of Sampled Signal',
            xaxis={'title': 'Sampled Normalized Frequency (f/fs)', 'showgrid': True, 'gridcolor': 'LightGrey'},
            yaxis={'title': 'Amplitude', 'showgrid': True, 'gridcolor': 'LightGrey'},
            margin={'l': 40, 'b': 40, 't': 40, 'r': 40},
            hovermode='closest',
            paper_bgcolor='white',
            plot_bgcolor='white',
            template='plotly_white'
        )
    }

    # Return this additional figure in the callback
    return sampled_signal_figure, sampled_fourier_transform_figure

    
    

if __name__ == '__main__':
    app4.run_server(debug=True, port=8051)

Φασματική Ανάλυση#

Για τον υπολογισμό της ενέργειας ή ισχύος της κυματομορφής \(x(t)\), ανάλογα με την περίπτωση σήματος, ισχύει

1.PNG

1.PNG

όπου για σήματα ισχύος \(S_Χ(f)\) είναι η πυκνότητα φάσματος ισχύος (Power Spectral Density – PSD) της \(x(t)\). Για σήματα διακριτού χρόνου που προκύπτουν από δειγματοληψία της \(x(t)\) με περίοδο \(T_s\), οι αντίστοιχες σχέσεις υπολογισμό της ενέργειας ή ισχύος γίνονται

1.PNG

1.PNG

Ένας απλός τρόπος να εκτιμηθεί η πυκνότητα φάσματος ισχύος της κυματομορφής \(x(t)\) είναι να ληφθεί ο DTFT των δειγμάτων του σήματος και μετά να υψωθεί στο τετράγωνο το μέτρο του αποτελέσματος. Αυτός ο εκτιμητής αποκαλείται περιοδόγραμμα (periodogram). Το περιοδόγραμμα ενός πεπερασμένου μήκους \(L\) σήματος \(x[n]\) ορίζεται ως

1.PNG

όπου \(X_d(\phi)\) o DTFT του σήματος. Με το μήκος \(L\) να τείνει στο άπειρο, το περιοδόγραμμα \(P_{xx}(f)\) τείνει στην πυκνότητα φάσματος ισχύος \(S_Χ(f)\). Ο υπολογισμός του περιοδογράμματος σε πεπερασμένο πλήθος συχνοτήτων \(kf_s/Ν\), \(k=0\), \(1\), \(…\) , \(Ν\) δίνει

1.PNG

όπου \(X_k\) και ο DFT της πεπερασμένου μήκους \(L\) σειράς δειγμάτων του σήματος. Η ισχύς του σήματος είναι τότε

1.PNG

όπου η τελευταία ισότητα προκύπτει από το θεώρημα Parseval, που για την περίπτωση του DFT εκφράζεται ως:

1.PNG

Στην ειδική περίπτωση περιοδικών σημάτων έχουμε

1.PNG

1.PNG

όπου \(X[k]\) οι συντελεστές του αναπτύγματος σε σειρά Fourier και \(T_o\) η περίοδος του σήματος.

Μέρος 3: Εφαρμογή Α#

# Part 1: Create the signal
Fs = 1000                        # Sampling frequency 1000 Hz
Ts = 1 / Fs                      # Sampling period
L = 1000                         # Length of signal (number of samples)
T = L * Ts                       # Duration of signal
t = np.arange(0, (L - 1) * Ts, Ts)  # Time vector

x = np.sin(2 * np.pi * 30 * t) \
    + 0.8 * np.sin(2 * np.pi * 80 * (t - 2)) \
    + np.sin(2 * np.pi * 60 * t)  # 60 Hz component

# Time domain plot of x
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=x, mode='lines', line=dict(color='#00CC96')))
fig.update_layout(title='Time domain plot of x', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='t (sec)',
                  yaxis_title='Amplitude',
                  template='plotly_white',
                  xaxis=dict(range=[0, 0.2]),
                  yaxis=dict(range=[-2.5, 2.5]))
fig.show()


# Fourier transform
def nextpow2(i):
    n = 1
    while n < i:
        n *= 2
    return n

N = nextpow2(L)                 # Length of Fourier transform
Fo = Fs / N                     # Frequency resolution
f = np.arange(0, N) * Fo        # Frequency vector
X = np.fft.fft(x, N)            # Compute DFT for N points

# Frequency domain plot of x (positive frequencies)
fig = go.Figure()
fig.add_trace(go.Scatter(x=f[1:N], y=abs(X[1:N]), mode='lines', line=dict(color='#1F77B4')))
fig.update_layout(title='Frequency domain plot of x', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='f (Hz)',
                  yaxis_title='Amplitude',
                  template='plotly_white')
fig.show()

# Shift frequencies to center
f = f - Fs / 2
X = np.fft.fftshift(X)

# Two-sided spectrum of x
f_shifted = f - Fs/2
fig = go.Figure()
fig.add_trace(go.Scatter(x=f_shifted, y=abs(X), mode='lines', line=dict(color='#1F77B4')))
fig.update_layout(title='Two sided spectrum of x', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='f (Hz)',
                  yaxis_title='Amplitude',
                  template='plotly_white')
fig.show()

# Calculate power
power = np.multiply(X, np.conj(X)) / N / L

# Periodogram
fig = go.Figure()
fig.add_trace(go.Scatter(x=f_shifted, y=power.real, mode='lines', line=dict(color='#1F77B4')))  # Use .real here if necessary
fig.update_layout(title='Periodogram', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='Frequency (Hz)',
                  yaxis_title='Power',
                  template='plotly_white')
fig.show()
# Part 2 Προσθέστε θόρυβο στο σήμα

# Συμπληρώστε τον κώδικα για τη δημιουργία του σήματος θορύβου n με τη βοήθεια της συνάρτησης randn.
# Το διάνυσμα θορύβου n θα πρέπει να είναι του ίδιου μεγέθους με αυτό της ημιτονοειδούς κυματομορφής x του πρώτου μέρους.
# Σχεδιάστε το σήμα θορύβου στο διάστημα από 0 έως 0.2 sec και κλίμακα σε από -2 έως 2.
# Υπολογίστε το περιοδόγραμμα του n και σχεδιάστε την πυκνότητα φάσματος ισχύος του σήματος θορύβου.
# Προσθέστε το σήμα θορύβου και το x για να λάβετε το σήμα με θόρυβο s.
# Σχεδιάσατε το σήμα με θόρυβο s στο πεδίο του χρόνου στην περιοχή 0 έως 0.2 sec 
# και κλίμακα από -2 έως 2 καθώς και το αμφίπλευρο φάσμα του.

rand_n = np.random.randn(np.size(x))

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=rand_n, mode='lines', line=dict(color='#00CC96')))
fig.update_layout(title='Time domain plot of n', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='t (sec)',
                  xaxis_range=[0, 0.2],
                  yaxis_title='Amplitude',
                  yaxis_range=[-2, 2],
                  template='plotly_white')
fig.show()

# Correction for N calculation using bitwise operator
N = 2^nextpow2(L)
Fo = Fs/N
f = (np.arange(0, N)) * Fo
f_shifted = f - Fs/2
rand_N = np.fft.fft(rand_n, N)
rand_N = np.fft.fftshift(rand_N)
power_n = np.multiply(rand_N, np.conj(rand_N)) / N / L

fig = go.Figure()
fig.add_trace(go.Scatter(x=f_shifted, y=power_n.real, mode='lines', line=dict(color='#1F77B4')))
fig.update_layout(title='Periodogram of n', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='Frequency (Hz)',
                  yaxis_title='Power',
                  template='plotly_white')
fig.show()

s = x + rand_n

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=s, mode='lines', line=dict(color='#00CC96')))
fig.update_layout(title='Time domain plot of s', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='t (sec)',
                  xaxis_range=[0, 0.2],
                  yaxis_title='Amplitude',
                  yaxis_range=[-2, 2],
                  template='plotly_white')
fig.show()

S = np.fft.fft(s, N)
S = np.fft.fftshift(S)

fig = go.Figure()
fig.add_trace(go.Scatter(x=f_shifted, y=np.abs(S), mode='lines', line=dict(color='#1F77B4')))
fig.update_layout(title='Two sided spectrum of s', title_x=0.5, title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='f (Hz)',
                  yaxis_title='Amplitude',
                  template='plotly_white')
fig.show()
# Part 3. Πολλαπλασιασμός σημάτων

# Συμπληρώστε τον κώδικα δημιουργίας ενός ημιτονοειδούς σήματος συχνότητας
# 100 Hz και πολλαπλασιάστε με το προηγούμενο σήμα s.
# Τα δύο σήματα θα πρέπει να είναι του ίδιου μεγέθους.
# Σχεδιάστε το αποτέλεσμα στο πεδίο του χρόνου στην περιοχή 0 έως 0.2 sec
# και κλίμακα από -2 έως 2 καθώς και στο πεδίο της συχνότητας
# χρησιμοποιώντας τη συνάρτηση fftshift.

z=np.sin(2*np.pi*100*t)
y= np.multiply(z,s)

# Time domain plot of y
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=y, mode='lines', line=dict(color='#00CC96')))
fig.update_layout(title='Time domain plot of y', title_x=0.5, 
                  title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='t (sec)', xaxis_range=[0, 0.2],
                  yaxis_title='Amplitude', yaxis_range=[-2, 2],
                  template='plotly_white')
fig.show()
           
N = 2^nextpow2(L)
Fo = Fs/N
f = (np.arange(0, N)) * Fo - Fs/2
Y = np.fft.fft(y, N)
Y = np.fft.fftshift(Y)

# Frequency domain plot of y
fig = go.Figure()
fig.add_trace(go.Scatter(x=f, y=np.abs(Y), mode='lines', line=dict(color='#1F77B4')))
fig.update_layout(title='Frequency domain plot of y', title_x=0.5, 
                  title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='f (Hz)', yaxis_title='Amplitude',
                  template='plotly_white')
fig.show()

Μέρος 4: Εφαρμογή Β#

Να γραφεί σε Python συνάρτηση φασματικής ανάλυσης, παρόμοια με την signal.welch(): θα δέχεται ως είσοδο διάνυσμα πραγματικού σήματος καθώς και τη συχνότητα δειγματοληψίας, \(F_s\), και θα σχεδιάζει τη μονόπλευρη φασματική πυκνότητα του σήματος στην περιοχή \([0-F_s/2)\). Το σήμα θα τεμαχίζεται σε τμήματα μήκους ίσου με τη δύναμη του \(2\) την πλησιέστερη στο \(1/8\) του συνολικού του μήκους, αλλά όχι μικρότερου από 256. Τα τμήματα θα είναι επικαλυπτόμενα κατά \(50\%\). Το τελευταίο τμήμα, εάν υπολείπεται σε μήκος των άλλων, θα αγνοείται. Θα υπολογίζεται με FFT το φάσμα κάθε τμήματος και θα λαμβάνεται η μέση τιμή όλων των τμημάτων. Η συνάρτηση να δοκιμαστεί με το σήμα του παραδείγματος 1.1 και να συγκριθεί το αποτέλεσμα με το αντίστοιχο της signal.welch().

def pwelch(x,Fs):                    
    Ts=1/Fs                    
    L=np.size(x)+1                 
    T=L*Ts                     
    N = 2^nextpow2(L)
    Fo=Fs/N                   
    f=np.arange(0,N)*Fo       
     
    window_size = nextpow2(np.size(x)/8)
    if (window_size<256):
        window_size=256
    windows = np.size(x)//(window_size//2)-1
    indexer = np.arange(window_size)[None, :] + (window_size//2)*np.arange(windows)[:, None]
    windowed_x = x[indexer]

    avg_pwr=0
    for window in windowed_x:
        window = window * np.hanning(np.size(window))
        L=np.size(window)+1                 
        T=L*Ts                     
        N = 2^nextpow2(L)
        Fo=Fs/N                   
        f=np.arange(0,N)*Fo
        window_fft=np.fft.fft(window,N)
        power=np.multiply(window_fft,np.conj(window_fft))/N/L
        avg_pwr=avg_pwr+power
    avg_pwr=avg_pwr/windows

    fig = go.Figure()
    fig.add_trace(go.Scatter(x=f[:N//2], y=avg_pwr[:N//2].real, mode='lines', line=dict(color='#1F77B4')))
    fig.update_layout(title='Periodogram pwelch()', title_x=0.5, 
                      title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                      xaxis_title='Frequency (Hz)', yaxis_title='Power',
                      template='plotly_white')
    fig.show()
    
    return f[np.arange(0,N//2)], avg_pwr[np.arange(0,N//2)]
Fs=500
f1,Pxx1 = pwelch(x,Fs)
f2,Pxx2 = signal.welch(x,fs=Fs)

fig = go.Figure()
fig.add_trace(go.Scatter(x=f2, y=Pxx2, mode='lines', line=dict(color='#1F77B4')))
fig.update_layout(title='Periodogram signal.welch()', title_x=0.5, 
                  title_font=dict(size=20, color='black', family="Arial, sans-serif"),
                  xaxis_title='Frequency (Hz)', yaxis_title='Power',
                  template='plotly_white')
fig.show()